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Fast Orthogonal transforms for pricing derivatives 

with quasi-Monte Carlo 

Christian Irrgeher* and Gunther Leobacher^ 


Abstract 

There are a number of situations where, when computing prices of financial 
derivatives using quasi-Monte Carlo (QMC), it turns out to be beneficial to apply 
an orthogonal transform to the standard normal input variables. Sometimes those 
transforms can be computed in time 0(nlog(n)) for problems depending on n input 
variables. Among those are classical methods like the Brownian bridge construction 
and principal component analysis (PCA) construction for Brownian paths. 

Building on preliminary work by Imai & Tan [3] as well as Wang & Sloan m, 
where the authors try to find optimal orthogonal transform for given problems, 
we present how those transforms can be approximated by others that are fast to 
compute. 

We further present a new regression-based method for finding a Householder 
reflection which turns out to be very efficient for a wide range of problems. We 
apply these methods to several very high-dimensional examples from hnance. 


1 Introduction 

Many simulation problems in finance and other applied fields can be written in the form 
E(/(V)), where / is a measurable function on M” and V is a standard normal vector, 
that is, X = (Vi,... ,Xn) is jointly normal with E(Vj) = 0 and E(VjVfc) = 6jk- It is a 
trivial observation of surprisingly big consequences that 

E{f{X))=E{f{UX)) (1) 

for every orthogonal transform U : M"" —?■ M”. While this reformulation does not change 
the simulation problem from the probabilistic point of view, it does sometimes make a 
big difference when quasi-Monte Carlo simulation is applied to generate the realizations 
of V. 

Examples are supplied by the well-known Brownian bridge and PCA constructions 
of Brownian paths which will be detailed in the following paragraphs. Assume that one 
wants to know E{g{B)) where i? is a Brownian motion with index set [0,T]. In most 
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applications this can be reasonably approximated by E,{g{BT,...,BTn)), where 5 ^ is a 

n n 

function on the set of discrete Brownian paths. 

There are three classical methods for sampling from {Bt, ..., Bnr) given a standard 

n n 

normal vector X, namely the forward method, the Brownian bridge construction and 
the principal component analysis construction (PCA). All of these constructions may be 
written in the form {Bt, ..., B^) = AX, where A is an n x n real matrix with 


AA = S := — min(j, k) 
\n 


n 

j,k=i n 


/ 111 . 
12 2 . 
12 3. 


1 2 3 


1 \ 
2 
3 

n 


For example, the matrix A corresponding to the forward method is 


A 


( 1 0 



V 1 1 


0 \ 
0 

1 / 


( 2 ) 


while PCA corresponds to A = VD, where S = VD‘^V~^ is the singular value decompo¬ 
sition of S. A corresponding decomposition for the Brownian bridge algorithm is given, 
for example, by Larcher, Leobacher & Scheicher [S]. 

It has been observed by Papageorgiou[TT] that AA^ = S if and only if A = SU for 
some orthogonal matrix U, so that every linear construction of {Bt, ..., B^) corresponds 

n n 

to an orthogonal transform of M”. In that sense the forward method corresponds to the 
identity, PCA corresponds to S~^VD and Brownian bridge corresponds to the inverse 
Haar transform, see [H]. 

Thus our original simulation problem can be written, as 


E(^(i?z, ...,Bt^))= ng^sx)) = E(/(X)) 

n n 

with f = g o S, and we interpret this as using the forward method. Consequently, the 
same problem using the Brownian bridge takes on the form K{f{H~^X)), where H is the 
matrix of the inverse Haar transform, and has the form K{f{S~^VDX)), with S, V, D 
as above, when PCA is used. 

As an application one can generalize the classical constructions of discrete Brownian 
paths to discrete Levy paths. See 00111 . 


There are some theories as to why an orthogonal transform might have the effect to 
make the problem more suitable for QMC. Caflisch et ah pQ introduce the concept of 
effective dimension of a function: consider a function g : M"" —M with finite variance 
w.r.t. normal distribution, that is E( 5 f(X)^) < cx) where X = {Xi,... ,Xn) is a vector of 
independent standard normal random variables. Then g may be written uniquely as the 
sum of functions g^ '■ E" E, u C {1 ,... ,77.}, where gu depends on the i-th coordinate 
only if 7 G M and where K{gu{X)) = 0 for all n 7 ^ 0 and K{gu{X)gy{X)) = 0 ior u ^ v, 
using the so-called ANOVA decomposition of g. Furthermore it holds 

n9(x)) = y v(9„(x)), 
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The effective dimension in the truncation sense at level a G (0,1) is then the smallest 
integer k such that 


V(j(V))(l-a)< V(j„(V)), 

0/«C{l. b) 

see [T]. Typically a is chosen as 0.01. Therefore, a function with effective dimension k 
is one that, in this sense, almost exclusively depends on the first k variables and which 
therefore is more suitable for QMC. This is confirmed by empirical evidence. Building on 
the concept of effective dimension of a function, Owen [njl gives definitions of effective 
dimensions of function spaces, thus connecting the concepts of effective dimension with 
that of tractability. 

Now one can turn this around and try to put as much variance as possible to the first 
few coordinates, by concatenating g with a suitable orthogonal transform. This is what 
has been done by Imai & Tan [3] and what we will do here, using a different approach. 
We shall see in Section 0] that empirical evidence also supports the conjectured efficiency 
of our method. 

However, there is also a disadvantage of that approach: the computation of the or¬ 
thogonal transform has a cost, which is in general of the order 0(?7,^). For large n this 
cost is likely to swallow the potential gains from the transform. We therefore concentrate 
on orthogonal transforms which have cost of the order 0(nlog(n)). 

Examples include discrete sine and cosine transform, Walsh and Haar transform as 
well as the orthogonal matrix corresponding to the PCA, see mM- 

Imai & Tan [3] propose an algorithm to find a good orthogonal transform in the sense 
that it puts as much variance as possible to the first few dimensions. They propose to 
take the first order Taylor expansion at some point X, i.e. 

g{X) » g(X) + _ x.). 

Then the contribution of the i-th component of X to Y{g{X)) is given by 

The columns of the orthogonal transform are chosen by solving optimization problems of 

the form 

fdg(AX) Y 

with I \A.i \I = 1 and = 0, j = 1,..., i — 1 

with Xi = {Xl ,..., X*, 0,..., 0). They suggest to perform this optimization only for the 
first few columns of the matrix A. In this paper we improve on their algorithm in various 
directions. In particular we find a good orthogonal transform that is fast in that it can 
be computed even in linear time. 

The remainder of the paper is organized as follows. Section [2] reviews basic proper¬ 
ties of Householder reflections and shows how they can be used to find fast versions of 
orthogonal transforms which put most variance on the first k variables. The main part 
of our article. Section [31 describes algorithms for finding fast orthogonal transforms using 
again Householder reflections. In contrast to the method of Imai & Tan [5] we do not 
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rely on differentiability. This makes the algorithm useful for barrier-type options. We 
further provide some theoretical results which indicate why the method serves to reduce 
the effective dimension. 

Section H] gives some numerical examples where the methods described earlier are 
applied to examples from hnance. We will see that the new methods described in Section [3] 
are among the best, both with regard to speed and accuracy. 

We provide an appendix where we compute certain expectations depending on the 
maximum of a Brownian path. This is useful for some of the numerical examples. 


2 Householder reflections 


We recall the dehnition and basic properties of Householder reflections from Golub & van 
Loan [2]. 

Definition 2.1. A matrix of the form 


U 


1-2 


vv 

v'^v 


where v G M”, is called a Householder reflection. The vector v is called the defining 
Householder vector. 


In the following proposition, ci denotes the hrst canonical basis vector in M”, ei = 

( 1 , 0 ,..., 0 ). 

Proposition 2.2. A Householder reflection have the following properties: 

1. If X ^ M”' is a vector then Ux is the reflection of x in the hyperplane spanjn}-*-. In 
particular, U is orthogonal and symmetric, i.e. U~^ = U. 

2. Given any vector a G M” we can find v & MA such that for the corresponding 
Householder reflection U we have Ua = ||a||ei. The computation of the Householder 
vector uses 3n floating point operations. 

3. The computation of the matrix-vector multiplication Ux uses at most An floating 
point operations. 

Proof. See Chapter 5.1 of Golub & van Loan 0 . □ 

Our main application of Householder reflections is the following: suppose we know 
that for a given integration problem E(/(X)) some orthogonal transform U reduces the 
effective dimension in the truncation sense to k, that is, almost all of the variance of 
fijjX) is captured by Xi,..., Xk, k <^n. 

Let U = (ui,... ,Un), that is, Uj is the j-th column of U. We can hnd Householder 
reflections Ui,... ,Uk such that Ui... UkCi = ui, £ = 1,..., k as follows: 

• Let f/i be a Householder reflection that maps ei to Ui. Ui also maps Ui to Ci. Since 
the vectors Uj are orthogonal we have eJ{UiU 2 ) = {UiUi)~^{UiU 2 ) = UiU 2 = 0 . 

• Therefore there exists a Householder reflection U 2 operating on the last n — 1 coor¬ 
dinates which maps 62 to U 1 U 2 . Thus UiU 2 ei = U^ei = Mi, UflJ 2 e 2 = UflJiU 2 = U 2 - 
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• Suppose Householder reflections Ui,... ,Uj have been constructed such that Ui.. .UjCi 
Ui, i= I,...,j. 

• Then there exists a Householder reflection operating on the last n — j coordi¬ 
nates which maps e^+i to Uj ... Uiitj+i. Then Ui... Uj+iei = iti, i = 1,..., j + 1. 

Write U = Ui.. .U^. By construction the first k columns of U coincide with those of 
IJ. Since, by assumption, Xi,... ,Xk capture almost all of the variance of f{UX), the 
same is true for f{UX). But for small k the computational cost for computing UX is of 
the order nk, as compared to general matrix-vector multiplication which occurs a cost of 
order 0(?7,^). 

Imai & Tan [3] and Wang & Sloan [T3] give examples for which they And good orthog¬ 
onal transforms U that reduce the effective dimension. However they do not specify how 
those transforms are applied. We propose to approximate them using the above method. 

However, the main topic of this paper is to present transforms that use only one 
Householder reflection. This will by detailed in the next section. 

3 Regression algorithm 

Let / : M” —> M be a measurable function with E(/(W)^) < cx) for a standard normal 
vector X. 

We want to approximate / by a linear function: 

f{x) ^ a^x + b 

where a G M"" and 6 G M. This can be done in different ways. For example, Imai & Tan 
13] take the first order Taylor expansion of /. 

In contrast, we take a “linear regression” approach, i.e. we minimize 

E ((/(W) - a^X - bf) min . (3) 

First order conditions give 

a,=E{f{X)X,), j = l,...,n and b = E{f{X)). 

Therefore, ([3]) minimizes the variance of the difference between / and the linear approxi¬ 
mation. So 

Y{f{X))=E{{f{X)-br) 

= E [{a^Xy + {f{X) - b - a^Xf) 

= I|a|r+V(/(X) -a^x) . 

That is, ||a||yV(/(X)) measures the proportion of variance captured by the linear ap¬ 
proximation. Now there exists a unique Householder reflection U that maps ei to a/||a||. 
With this transform we have aYUX = ||a||e7X = ||a||Xi and therefore 

E(/(X)) = EUiUX)) = E (a^UX + {f{UX) - a^UX)) 

= E(||a||Xi + (/(f/X)-||a||W)) • 
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Therefore the linear part of the integration problem depends on the parameter Xi alone. 
Now, if the linear part constitutes a large part of the integration problem then we have 
succeeded in putting a large fraction of the variance into the hrst coordinate by composing 
/ with U. 

Algorithm 3.1. Let be independent standard normal variables. Let f be a 

function f : M” —)• M. 

1. Oj := E{Xjf{X)) for j = 1,... ,n; 

2. if ||a|| = 0 define U = I and go to f.; 

3. else let U be a Householder reflection that maps Ci to a/||a|l; 


4- Compute E{f{UX)) using QMC. 


A drawback of the algorithm is that in general the computation of the expectations 
in step 1 is no easier than the original problem. In some cases the expectation can be 
computed explicitly, though usually in that case also the original problem has an explicit 
solution. 

Example 3.2. f{X) = Yfk=i'^k^'^V {Ylf=i{ck,jXj + dfcj)). It is easily verihed that, with 0 
denoting the standard normal density, 0(a;) = exp(— 


/ exp(cx + d)(j){x)dx = exp fc^/2 + d^ , 

'J M 

/ X exp(ca; + d)(j){x)dx = cexp (cf /2 + d 

^ M. 


Therefore it holds that 

/ CO 

f{x)Xi(j){Xi) . . . (j){Xn)dXi ...dXr. 

-CO 


' —CO 

m 


= ^k,iWk exp ^ (cfcj/2 + dkj 
k=l \j=i 


Let us hnd out how much of the variance of f{UX) is captured by ||a||Xi: 

We write Wk := tCfcexp(X]"=i(c|j/2 + dkj)). Then 

n / m 

||a||^ = XI H '^kCk^i 

i=l \k=l 
mm n 

^ Wk^Wk2 ^ ^ Ck^^iCk2,i 

fci=l k2=l *=1 

m m 

'y ' y ' 'djkiUJk2Cki,k2 7 (4) 

fci=l k2=l 

where Ck,,k 2 ■= Ei=i Cki,iCk 2 ,i- 

On the other hand, it is easy to see that E(f(X)) = ET=i'^k and E(/(A)^) = 
Wkj^Wk 2 e^’‘^’^^. Therefore we get for the variance of f{UX) 

Y{f{UX)) = V(/(A)) = E(/(A)^) - E(/(A))2 

m m 

= X X Wk.Wk^ie'''^^’^^ - 1). (5) 

ki=l k2=l 
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Let us try some special values that are related to Asian options: 


1 

TTi Tl, 

n 


= a 








with r, cr, T > 0, At = ^. For this choice we get Wk = and Cfc^.fca = . 

For large n the sums in equations (jl]) and ([S]) can be approximated by corresponding 
integrals such that 


||a||^ cr^T [ [ mm{x,y)dxdy 

Jo Jo 

2 4e’'^ + 2e^^^rT - + 1) 

“ ^ 2r3T2 

V(/(X))2^ -l)dxdy 

Jo Jo 

2e’’^ (2rcr^ + a^) + (2r^ + 3rcr^ + a^) + ra^ + 

r‘^T‘^[r + cr^)(2r + a^) 


Table [T] shows the fraction 


v(/A))-lbll 

v(/A)) 


for a few values of r, a and T 


1 . 


r\a'^ 

0.01 

0.02 

0.03 

0.04 

0.1 

0.0025 

0.0051 

0.0076 

0.0101 

0.2 

0.0026 

0.0051 

0.0077 

0.0103 

0.3 

0.0026 

0.0052 

0.0078 

0.0104 


Table 1: for T = 1 and different values for r, 


It can be concluded that in this example almost all of the variance of f{UX) is captured 
by Xi. □ 

In general we cannot expect that E(/(X)Xj) can be compnted explicitly. Of course 
it is an option to compnte E(/(X)Xj) using (quasi-)Monte Carlo, though it is unlikely 
that this will lead to small overall computing times. But quite frequently, especially in 
hnancial applications, a problem can be written in the form, f{X) = g{h{X)), where 
E(h(X)Xj) can be compnted and h is some relatively simple fnnction h : M — > M. 

Algorithm 3.3. Let Xi,..., Xn be independent standard normal variables. Let f be a 
function f : R"' —> R, which is of the form f = goh where h : R” —)■ R and g : R —)■ R. 

1. Oj ;= E(Xjh(X)) for j = 1 ,..., n; 

2. if ||a|| = 0 define U = I and go to f.; 

3. else let U be a Householder reflection that maps ei to a/||a|l; 

4- Compute E(/(t/X)) using QMC. 

Without additional assumptions on the functions h and g there is no guarantee that U 
gives better convergence. Nevertheless there are practical examples where this algorithm 
gives excellent results. 
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Example 3.4. Consider an arithmetic average value option written on some underlying S, 

f{X) = e-^^max - iC, 0^ , 

and 

S^AX) = Soexp " y) • 

Here we have f{X) = g{h{X)), where g{s) = max(s — K, 0) and h is like in Example 
Owith m = n,Wk = ^Sq, Ckj = dkg = - ^)lj<k ■ □ 

Write Y := ||a||Xi, Z := h{UX) — ||a||Xi. Then Y,Z are uncorrelated, 

E{YZ) =E(h(17X)||a||e7X) - ||af = E(h(t/X)||a||(t/ei)^17X) - ||a||^ 
=E{h{UX)a'^UX) - ||af = a'^E{h{UX)UX) - Hajl^ 

=a^E(h(X)X) - ||a|p = a^a - ||a|p = 0 . 

Further, E(F) = 0, such that E(F)E(Z) = 0 as well, and therefore Cov(y, Z) = 0. 

Theorem 3.5. Let f,g,h,U,Xi,...,Xn be like in Algorithm VJ.IA Write again Y : = 
||a||Xi, Z := h{UX) - ||a||Xi. 

Then Y{f{UX)) = V(E(^ 7 (F + Z)\Y)) + V(r?(F + Z) - E{g{Y + Z)\Y)). 

Proof. We write Y =_E{g{Y + Z)|F)_and Z = g{Y + Z) - E{g{Y + Z)\Y), so that we 
have to show V{Y + Z) = \{Y) + V(Z). To that end it is sufficient to prove that Y and 
Z are uncorrelated: 

EiYZ) =E(E(r?(r + Z)\Y)g{Y + Z)) - E(E(^ 7 (r + Z)\Y)E{g{Y + Z)\Y)) 
=E(E(E(r/(F + Z)\Y)g{Y + Z)\Y)) - E(E(^ 7 (r + Z)\Yf) 

=E(E(( 7 (F + Z)\Y)Eig{Y + Z)\Y)) - E(E(( 7 (y + ZfY^) = 0 . 

Since E{Z) = 0, we have E{Y)E{Z) = 0 = E(YZ). □ 

We consider a special case that will rarely occur in practice but which gives a flavor of 
the best result possible. Assume that g is Lipschitz continuous with constant L. Suppose 
further that Y and Z are not only uncorrelated, but even independent. 

Denote by Fy, Fz the cumulative probability distribution functions of Y and Z, re¬ 
spectively. Using independence we get 

E{g{Y + Z)\Y)= [ g{Y + OdFziC). 

Jm. 

Noting that E{g{Y + Z) — E{g(Y + Z)\Y)) = 0 we thus get 

V(g(r + Z) - E(g(r + Z)\Y)) = E (( 9 (r + Z) - E( 9 (r + Z)\Y)f ) 

= E ((£ {g{Y + Z) - 9 (r + OyFziof] 
<1^(1 {g(Y + z) - 9(r + C))" dFzio'j 
<e(^L^ f(Z-(fdFz(0'j 
< L^E - 2ZE(Z) + E{ZA) = 2L‘^N{Z ), 




where we also have used the Cauchy-Schwarz inequality. Thus with Theorem IT5I we get 

W{f{UX)) - Y(E{g{Y + Z)\Y)) < 2L‘^Y{Z) 

that is, 

V(/(f/X)) - V(E(/(f/X)|Xi)) < 2L^{Y{h{UX)) - ||a|n. 

So in this situation, if Xi captures a large fraction of the variance of h{UX), then Xi also 
captures a large fraction of the variance of f{UX) provided that the Lipschitz constant 
L is not too big. 

We can also think of a variant of Algorithm [T3] for slightly more complicated functions. 
We have been inspired by Wang & Sloan na, where the authors consider functions of the 
form f{X) = g{wJX,... ,w^X) and show, that there is an orthogonal transform that 
makes this function m-dimensional. We give a slightly modihed version of their argument 
which guarantees that the orthogonal transform is also fast to compute for small m, that 
is for m < log(n). 

Let f{X) = g{wJX ,..., w^X) for tci,..., Wm G M”. We may assume that wi is not 
the zero vector. Let f/i be a Householder reflection which maps (1, 0,..., 0) to Then 
wJUiX = ||tci||(l, 0 ,..., 0 )'''X = ||tci||Wi and therefore 

f{U,X) = g{\\wi\\Xi, {U,W2VX ,..., (UiwJ^X ). 

Next we write {UiWkyX = {UiWk)JXi + {UiWk) 2 ...n^ 2 ...n- That is, 

f{U,X) = g{X^, wJX2...n, wlX2...n) 

for some W 2 , ■ ■ ■, Wm G Assuming that ti )2 7 ^ 0, let U 2 be the Householder reflection 

of that maps ( 1 , 0 ,..., 0 ) to tC 2 /||'u) 2 || and let 


Then U 2 is a Householder reflection of R” and 

/(Hif/aX) = |(Ai, A2, ..., . 

for some ... ,Wn G R”“^. Proceeding that way one arrives at 

f{Ur---U^X) = g{X^,X2,...,X^) 

for some rh <m (We may have m < m if at some stage all remaining Wk are zero). 

We propose a similar procedure for an integration problem of the form f{X) = 
g{hi{X), h 2 {X),..., hm{X)) where K{hj{X)Xk) can be computed explicitly (or at least 
efficiently). 

Algorithm 3.6. Let Xi,..., Xn be independent standard normal variables. Let f be a 
function f : R” —> R, which is of the form f = g o h where h : R"" —> R™' and 
g : R™' —R. 
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1. Start with k = 1 and U = I; 

2. af := ¥.{Xjhk{UX)) forj = k,..., n; 

3. af'’ ■= 0 for j = 1,... ,k - 1; 

4- if = 0 define = I and go to\^ 

5. else let be a Householder reflection that maps Ck to 

6. U = UU^’^fl 
1. k = k + 1; 

8. while k <m, go back to\^ 

9. Compute ¥.{f{UX)) using QMC. 

We will give a numerical example in Section 01 


4 Numerical tests 

In this section we will apply our method to examples from mathematical finance. 

Asian option 

The first numerical example we give is the evaluation of an Asian call option with discrete 
arithmetic average in the Black-Scholes model, which has been discussed previously. Since 
the payoff function / is of the form goh with g and h as in Example 13.41 we apply Algorithm 
13.31 to the integration problem E{f{X)) where the vector a follows from Example 13.21 i.e. 
for every i = 1,... ,n 



For the quasi-Monte Carlo simulation we use a Sobol sequence of dimension n = 250 with 
a random shift and we have Sq = 100, K = 100, r = 0.04, a = 0.2 as well as T = 1. 
We compute the standard deviation based on 32 batches for N sample paths, where 
the number of sample paths ranges from 2^ to 2^^. Note that the standard deviation is 
different from the RQMC standard deviation defined in L’Ecuyer & Munger [Sj. 

In Figure[T]we compare the regression method with the forward method, the PCA con¬ 
struction and the LT method of Imai and Tan. We see that PCA, LT and the Regression 
method yield similar results, but all of the three outperform the forward method. Note 
that the regression method can be applied in 0{n). Thus we can achieve the efficiency 
of the PCA with the regression method with lower computational costs. Moreover, it is 
interesting that the LT method and regression method yield nearly the same results. 

The computation time required to price the Asian option using quasi-Monte Carlo 
integration with 2^^ paths is given in Table El Note that PCA is implemented using the 
discrete sine transform as discussed in Leobacher [S]. The LT method is implemented 
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Figure 1: Asian option (left) and Asian Basket option (right): Standard deviation of 32 
runs on a log 2 -scale 



Forward 

PCA 

LT 

Regression 

time (sec) 

0.08 

0.64 

1.94 

0.15 


Table 2: Computation times for pricing the Asian option 


such that only the hrst 25 columns are optimized and then the orthogonal transform is 
completed using Householder reflections as we suggested in Section [2J 

Furthermore it should be mentioned that the regression method as well as the LT 
method produce an overhead caused by determining the orthogonal transform. Neverthe¬ 
less the overhead time is rather small and is negligible for a large sample size. 

The computation times of the subsequent numerical examples are similar to the result 
regarding the Asian option. 


Asian basket option 

We consider an Asian basket call option with arithmetic average and a basket consisting 
of m assets, an example taken from Imai & Tan [3]. The i-th asset of the basket 
(i = 1 ,..., m) is given by 



Ff^ exp 



k- + 
n 


where is the current price of the i-th asset, r is the risk-free interest rate, ai is the 
volatility of the i-th asset and B = {B^^\ ..., B^'^'^) is an m-dimensional Brownian motion. 
The correlation between B^^'^ and B^^'> is denoted by pjk- The payoff function of the Asian 
basket option is given by 


f{X) = max 



ni lb > 

Y.T.stl(x) - 


2 = 1 k = l 


1 
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where 


( nm 2 rp 

J2C(k-i)„,+ijXj + {r - ^)k- 

1 Z Ti 

j=i 

and where C is an mn x mn—matrix with CC^ = S := i?® S and i? is an m x m—matrix 
with Rii = ^JT/naf for all i and Rij = ^Jt/ n pijaiaj for i ^ j. Note that the discussion 
of the previous sections also holds for a discrete Brownian path with covariance matrix 
S. Since the problem is of the form f{X) = g{h{X)), Algorithm Id.dl can be applied. 
Since the function h is of the form considered in Example Id.21 we can compute the 
corresponding vector a analytically. Furthermore, notice that the PCA construction can 
be computed in this example efficiently by using the orthogonal transform ViDi ® V 2 D 2 
where ViDlV^ = R and V 2 DlV:^ = S. 

The parameters are T = 1, r = 0.04, K = 100 and pjk = 0.05 for j 7 ^ k. Moreover, the 
volatility of the 10 assets is equally spaced from 0.1 to 0.3 and we assume that = 100 
for all i = 1,..., m. Since we simulate every asset at 250 time points, we take a Sobol 
sequence in dimension n = 2500 with a random shift. In Figure [T] we can observe the 
standard deviation based on 32 batches of the forward method, the PCA construction, 
the LT method and the regression method for N sample paths with N up to 2^^. 

Digital barrier option 

A digital (up-and-in) barrier option is a derivative which pays 1 if the underlying asset 
breaks through a barrier u on the time interval [0, T] and pays 0 otherwise. We intend to 
price the option in a discrete Black-Scholes model, where the path of the stock is given 
by S' = {Si, ...,Sn) with 

Sk{X) = So exp (J^- 

with current stock price Sq, interest rate r, volatility a, Brownian path B = {B^t)'^^^ 

_ n 

where B^,t = \/-Yl’j=iXj and standard normal vector X = {Xi,..., Xn). Hence, the 

n \ Tl J 

payoff function h of the digital barrier option is 

h{X) = 

1 max/j— 1 ,...,n ‘5'fc (^ 

which leads us to an integration problem of the form E(exp(—rT)h(i?)). We can use 
Algorithm [23] for solving this problem and therefore, we have to compute a* = E(h(X)Xj) 
for i = 1,... ,n. In the appendix we show how to calculate this expectation for a function 
depending on the maximum of a Brownian motion with drift u. We can adjust our problem 
by 

max Sh > u 

k=l,...,n 


max So exp ( (r —— ] k -h aB,T ] > u 

k=l,...,n \ \ 2 j n " / 


max 


+ B,r_ > MhA 


a 


n 


a 


max B'j^T > u 


k=l,...,n 
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with = ut + Bt, V = ^ and u = With ([7j) we get that the vector a in 

Algorithm 13.31 can be approximated by 

a ^ — v^T/n 7 1 

where S is given by ([2]), /3 = (/3i,..., f^nY with /?* = E(lmaxo<,<TBr>«^r^)) 7 = E(lmaxo<,<TB.">«) 

n 

and 1 = (1,..., 1)^. The computation of f3i with i = 1,... ,n can be reduced to a 1- 
dimensional integration problem using flTUD with / = id® and t = and formula ([H]) with 
/ = 1 simplihes 7 . Consequently, we end up with 1-dimensional integrals which can be 
evaluated efficiently with an adaptive quadrature rule. 

For the numerical test we use a Sobol sequence of dimension n = 2000 with a random 
shift and the parameter set is chosen as Sq = 100, u = 110, r = 0.04, a = 0.2 and T = 1. 

The number of sample paths N ranges from 2^ to 2^“^ and we compute the standard 
deviation for those N based on 32 batches. Since it is not clear how to apply the LT 
method of Imai and Tan to barrier options, we compare the regression method with the 
forward method and the PCA construction only. In Figure E] we can observe that the 
difference between the forward method and the PCA is smaller than in the previous 
examples. Furthermore, we see that the regression method is slightly behind the PCA, 
but this seems to be the best we can achieve by linear approximation. 
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- * - Regression 






8 10 12 14 
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10 12 14 


Figure 2: Digital barrier option (left) and Asian barrier option (right): Standard deviation 
of 32 runs on a log 2 -scale 


Asian barrier option 

The last example we provide is an Asian (up-and-in) barrier option by which we mean that 
the payoff of the option is similar to an Asian option as in the hrst numerical example, but 
is paid only if the underlying asset breaks through an upper barrier u. The corresponding 
function is then given by 

f{X) = exp(-rT)lmaxfe^i,...,„5fe(x)>nmax [^-S'fc(X) - K, Oj 

where Sk{X) is as in ([H]) for k = 1,... ,n. Since the function / is of the form f{X) = 
g{hi{X),h2{X)) with g{x,y) = exp(-rT)a;|/, hi{X) = lmaxfe^i,...,„Sfe(x:)>« and h2{X) = 
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max(X]fc=i ^-S'fc(X) — 0), we apply Algorithm ESI with m = 2 to the problem. The 

computation of the vectors and is already discussed in the examples above, i.e. 
is related to the digital barrier option and corresponds to the Asian option. 

The numerical test is based on 32 batches and we again compare the standard deviation 
of the forward method, the PCA construction and the regression method for various 
numbers of sample paths N, ranging from 2^ to 2^“^. Moreover, we use a Sobol sequence 
in dimension n = 1000 with a random shift and the parameters are Sq = 100, K = 
100, M = 110, r = 0.04, a = 0.2 and T = 1. Figure [2] shows that the regression method 
yields slightly better results than the PCA and that the forward method is behind the 
other two approaches. 


Appendix: Regression for the maximum 


We give the computations needed for examples of barrier type, that is we want to compute 
E(h(X)Xj) where h is some function of the maximum of a discrete Brownian path with 
drift z/, i.e. 


h{X) 


h max 

\ k 


B, 



5 


and where 

n 



Xj. We make the approximation 


E h max B^ + u 


kT' 


n 



E ( h ( max [Bg + us) 
V \o<s<B 


y/n ( i?iT — B{i- 


1)T 


= E (ft (m|« B;) ^/S (b| - - A)) . (7) 


where B'^ denotes Brownian motion with drift z/ G M, i.e. B)) := Bt + ut, t > 0. Moreover, 
let Mj^j, := m8iXt<s<T and := Mq ^. At hrst we compute K^lM^^ufiB)))) for given 
M > 0 and measurable / with E(|/(i?^^)|) < cx). Then we show how the expectation for 
more general h can be computed using the hrst result. 

We start with a simple calculation for a Brownian motion B with drift 0 and let 
Mt := A/°. For u > 0 we get, using the rehection principle for Brownian motion. 


^i)-Mt>uf{.Bt)) — ¥.{lMt>BXBt>uf\Bt)) +¥,{lMt>vXBt<uf{Bt)) 

= ¥.{lBt>uf{Bt)) + E(lsj>„/(2n — Bt )). 


Next we make a Girsanov-type change of measure such that under the new measure 
Q the Brownian motion B'^ with drift becomes a standard Brownian motion. So with 


dQ 

dP 




> that is g = e 




= + EQ(lBj^>„/(2n - 

= E(lBJ^>„/(5n) + ^Q{l-B~£>uf{2u + 

= E(lBr>n/(i?r)) + e2-E(lB,^<_„/(2n + 5^)) • (8) 
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denote the standard 


The next step is to consider E(1 m^>u/(-B^)) for t <T. Let {J^t}o<t<T 
hltration of B. 

= E,(^f {B’j')'E{lM^>u + 

= E[f{B’^)lMr>u) + E(lMr<uf{B’^)mM^^^^>u\J^t)) ■ (9) 

We have already computed the hrst term. For the second term we note that by the 
Markov property of Brownian motion, 




We can use our earlier result ([HD with f{x) = l to obtain 


E(1 




Rfo = $ 


B^ — u — p{T — t) 


VT^ 


(1 + . 


Let us write g{u,x) := Then, using (|HD and (|HD we obtain 


=E{f{B^)g{u, B^)) + E{f{B’^)lMr>u) - E{lMr>uf{B^)g{u, B^)) 

=E{f(BMu, Bn) + E (l„->„/(Br)(l - ff(u, Bn)) 

=E(f(Bng(n, Bn) +E (iB~>„f(Bn(i - »(«, sr))) 

+ e^“‘'E (lB.<_„/( 2 t. + Bn(l - 9 {n, 2 u + Bn)) ■ (10) 

Note that the expectations can be computed explicitly for suitable /. 

We can also use dTUp to compute E(h(M^)f(Bf)) for h differentiable and h(0) = 0 
and such that the expectations all converge absolutely: 


E(h(M^)f(Bn) =E(E(h(M^)lBnf(Bn) 

=E^“ fo(n)E(lM^>„|i?r)c^/ra) 

POO 

= h\u)E{E{lM!^>^\B'i)f{B'i))du 

POO 

= J^ h\u)E{lM)f>uf{B‘^))dii. 
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